

##### Load Data

	setwd("C:/Users/eas24f/OneDrive - Florida State University/Research/Inertia RESTAT/Data")
	#demand_data <- read.csv(file="ca_julia_data_AUG032019_small.csv",header=TRUE)
	households <- get(load("ca_household_data_AUG022019"))
	outside_logit <- get(load("sipp_logit"))
	gc()

	# Drop uninsured households
	households <- households[!is.na(households$plan_number_nocsr),]
	
	plan_data <- read.csv("ca_plan_data2.csv") # Covered California plan data by year and rating area
	plan_data$Plan_ID_Order <- as.character(plan_data$Plan_ID_Order)
	rownames(plan_data) <- paste(plan_data$Plan_Name,plan_data$ENROLLMENT_YEAR,
		plan_data$region,sep="")
	plan_data <- plan_data[plan_data$Plan_ID_Order,]
	plan_data <- plan_data[plan_data$ENROLLMENT_YEAR < 2019,]
	
	risk_score_reg <- get(load("rs_reg_full"))
	
##### Household AV variable

	plan_data$AV <- 0.6
	plan_data[plan_data$metal_level == "Silver","AV"] <- 0.7
	plan_data[plan_data$metal_level == "Silver - Enhanced 73","AV"] <- 0.73
	plan_data[plan_data$metal_level == "Silver - Enhanced 87","AV"] <- 0.87
	plan_data[plan_data$metal_level == "Silver - Enhanced 94","AV"] <- 0.94
	plan_data[plan_data$metal_level == "Gold","AV"] <- 0.8
	plan_data[plan_data$metal_level == "Platinum","AV"] <- 0.9
	households$AV <- plan_data[households$plan_id,"AV"]
	
##### Delete CSR plans and reassign households plan_ID to silver

	plan_data$unique_id <- paste(plan_data$Plan_Name,plan_data$region,plan_data$ENROLLMENT_YEAR,sep="")
	plan_data$unique_id_nocsr <- paste(plan_data$Plan_Name_NOCSR,plan_data$region,plan_data$ENROLLMENT_YEAR,sep="")
	households$unique_id_nocsr <- plan_data[households$plan_id,"unique_id_nocsr"] 
	plan_data <- plan_data[!duplicated(plan_data$unique_id_nocsr),]
	rownames(plan_data) <- plan_data$unique_id_nocsr
	
##### Create Plans PMT object

	plan_data$Insurer <- plan_data$Issuer_Name
	plan_data$Insurer <- as.character(plan_data$Insurer)
	plan_data[plan_data$Insurer == "Anthem Blue Cross","Insurer"] <- "Anthem"
	plan_data[plan_data$Insurer == "Blue Shield","Insurer"] <- "Blue_Shield"
	plan_data[plan_data$Insurer == "Chinese Community","Insurer"] <- "Chinese_Community"
	plan_data[plan_data$Insurer == "Contra Costa Health Plan","Insurer"] <- "Contra_Costa"
	plan_data[plan_data$Insurer == "Health Net","Insurer"] <- "Health_Net"
	plan_data[plan_data$Insurer == "LA Care","Insurer"] <- "LA_Care"
	plan_data[plan_data$Insurer == "Western Health","Insurer"] <- "Western_Health"
	
	plan_data$Year <- plan_data$ENROLLMENT_YEAR
	plan_data$market_year <- paste(plan_data$region,plan_data$Year,sep="")
	plan_data$Metal_Level <- plan_data$metal_level
	plan_data$Type <- plan_data$PLAN_NETWORK_TYPE
	plan_data$Plan_Name_Small_NOHSACAT <- as.character(plan_data$Plan_Name_Small_NOHSACAT)
	plan_data$insurer_small <- plan_data$Insurer
	plan_data[!plan_data$insurer_small %in% c("Anthem","Blue_Shield","Health_Net","Kaiser"),"insurer_small"] <- "Small_Insurer"
	plan_data$Anthem <- as.numeric(plan_data$insurer_small == "Anthem")
	plan_data$Blue_Shield <- as.numeric(plan_data$insurer_small == "Blue_Shield")
	plan_data$Health_Net <- as.numeric(plan_data$insurer_small == "Health_Net")
	plan_data$Kaiser <- as.numeric(plan_data$insurer_small == "Kaiser")
	
	#keep_fields <- c("pmt_name","Plan_Name_Small_NOHSACAT","ENROLLMENT_YEAR","region","Metal_Level","PLAN_NETWORK_TYPE","Premium",
	#	"insurer_small","Anthem","Blue_Shield","Health_Net","Kaiser")
	plan_data$pmt_name <- paste(plan_data$Plan_Name_Small_NOHSACAT,plan_data$region,plan_data$ENROLLMENT_YEAR,sep="")
	plan_data$pmt_premium <- by(plan_data$Premium,plan_data$pmt_name,mean)[plan_data$pmt_name]
	plan_data$HMO <- as.numeric(plan_data$PLAN_NETWORK_TYPE == "HMO")
	
	plans_pmt <- plan_data[!duplicated(plan_data$pmt_name),]
	rownames(plans_pmt) <- plans_pmt$pmt_name
	
	plans_pmt$plan_name <- plans_pmt$Plan_Name_Small_NOHSACAT
	plans_pmt$year <- plans_pmt$ENROLLMENT_YEAR
	plans_pmt$rating_area <- plans_pmt$region
	plans_pmt$market_year <- paste(plans_pmt$rating_area,plans_pmt$year,sep="")
	plans_pmt$HMO <- as.numeric(plans_pmt$PLAN_NETWORK_TYPE == "HMO")
	plans_pmt$Insurer <- plans_pmt$insurer_small
	
	#plans_pmt <- plans_pmt[,keep_fields]
	
	# Need to add risk score and market share
	
		# Input plan into household object
		households$plan_name <- plan_data[households$unique_id_nocsr,"Plan_Name_Small_NOHSACAT"]
		households$pmt_name <- paste(households$plan_name,households$rating_area,households$year,sep="")
		households$market_year <- paste(households$rating_area,households$year,sep="")
	
		# Identify plans with positive enrollment
		positive_plans <- plans_pmt[plans_pmt$pmt_name %in% households$pmt_name,"pmt_name"]
		pd_positive_plans <- unique(households$unique_id_nocsr)
		
		# Create market share variable
	
		fields <- c("demand","market_demand","market_share")
		plans_pmt[,fields] <- 0
		plans_pmt[positive_plans,"demand"] <- by(households$weight,households$pmt_name,sum)[positive_plans]
		plans_pmt$market_demand <- by(households$weight,households$market_year,sum)[plans_pmt$market_year]
		plans_pmt$market_share <- plans_pmt$demand/plans_pmt$market_demand
		
		plan_data[,fields] <- 0
		plan_data[pd_positive_plans,"demand"] <- by(households$weight,households$unique_id_nocsr,sum)[pd_positive_plans]
		plan_data[is.na(plan_data$demand),"demand"] <- 0
		plan_data$market_demand <- by(households$weight,households$market_year,sum)[plan_data$market_year]
		plan_data$market_share <- plan_data$demand/plan_data$market_demand
		
		# Create risk scores
		
		fields <- c("Intercept","Silver","Gold","Platinum","share_18to34","share_35to54","share_250to400","share_gt400",
			"share_male","share_family","share_Asian","share_Black","share_Hispanic","share_Other_Race")
		
		plans_pmt[,fields] <- 0
		plans_pmt[,"Silver"] <- as.numeric(plans_pmt[,"Metal_Level"] == "Silver")
		plans_pmt[,"Gold"] <- as.numeric(plans_pmt[,"Metal_Level"] == "Gold")
		plans_pmt[,"Platinum"] <- as.numeric(plans_pmt[,"Metal_Level"] == "Platinum")
		
		plan_data[,fields] <- 0
		plan_data[,"Silver"] <- as.numeric(plan_data[,"Metal_Level"] == "Silver")
		plan_data[,"Gold"] <- as.numeric(plan_data[,"Metal_Level"] == "Gold")
		plan_data[,"Platinum"] <- as.numeric(plan_data[,"Metal_Level"] == "Platinum")
		
		
		plans_pmt$AV <- 0.7
		plans_pmt[plans_pmt$Metal_Level %in% c("Bronze","Minimum Coverage"),"AV"] <- 0.6
		plans_pmt[plans_pmt$Metal_Level %in% c("Gold"),"AV"] <- 0.8
		plans_pmt[plans_pmt$Metal_Level %in% c("Platinum"),"AV"] <- 0.9
		
		plan_data$AV <- 0.7
		plan_data[plan_data$Metal_Level %in% c("Bronze","Minimum Coverage"),"AV"] <- 0.6
		plan_data[plan_data$Metal_Level %in% c("Gold"),"AV"] <- 0.8
		plan_data[plan_data$Metal_Level %in% c("Platinum"),"AV"] <- 0.9
		
		
		#plans_pmt[plans_pmt$Metal_Level %in% c("Minimum Coverage","Bronze"),"Metal_Level"] <- "Bronze"
		#plans_pmt[plans_pmt$Metal_Level == c("Bronze"),"AV_Demean"] <- 0.6 - 0.7
		#plans_pmt[plans_pmt$Metal_Level == c("Gold"),"AV_Demean"] <- 0.8 - 0.7
		#plans_pmt[plans_pmt$Metal_Level == c("Platinum"),"AV_Demean"] <- 0.9 - 0.7	
		
		
		plans_pmt$Intercept <- 1
		plans_pmt[positive_plans,"share_18to34"] <- by(households$weight * (households$perc_18to25 + households$perc_26to34),
			households$pmt_name,sum)[positive_plans]/plans_pmt[positive_plans,"demand"]
		plans_pmt[positive_plans,"share_35to54"] <- by(households$weight * (households$perc_35to44 + households$perc_45to54),
			households$pmt_name,sum)[positive_plans]/plans_pmt[positive_plans,"demand"]
		plans_pmt[positive_plans,"share_250to400"] <- by(households$weight * (households$FPL > 2.5 & households$FPL <= 4),
			households$pmt_name,sum)[positive_plans]/plans_pmt[positive_plans,"demand"]
		plans_pmt[positive_plans,"share_gt400"] <- by(households$weight * (households$FPL > 4),
			households$pmt_name,sum)[positive_plans]/plans_pmt[positive_plans,"demand"]
		plans_pmt[positive_plans,"share_male"] <- by(households$weight * households$perc_male,
			households$pmt_name,sum)[positive_plans]/plans_pmt[positive_plans,"demand"]
		plans_pmt[positive_plans,"share_family"] <- by(households$weight * (households$household_size > 1),
			households$pmt_name,sum)[positive_plans]/plans_pmt[positive_plans,"demand"]
		plans_pmt[positive_plans,"share_Asian"] <- by(households$weight * households$perc_asian,
			households$pmt_name,sum)[positive_plans]/plans_pmt[positive_plans,"demand"]
		plans_pmt[positive_plans,"share_Black"] <- by(households$weight * households$perc_black,
			households$pmt_name,sum)[positive_plans]/plans_pmt[positive_plans,"demand"]
		plans_pmt[positive_plans,"share_Hispanic"] <- by(households$weight * households$perc_hispanic,
			households$pmt_name,sum)[positive_plans]/plans_pmt[positive_plans,"demand"]
		plans_pmt[positive_plans,"share_Other_Race"] <- by(households$weight * households$perc_other,
			households$pmt_name,sum)[positive_plans]/plans_pmt[positive_plans,"demand"]
		
		
		plan_data$Intercept <- 1
		plan_data[pd_positive_plans,"share_18to34"] <- by(households$weight * (households$perc_18to25 + households$perc_26to34),
			households$unique_id_nocsr,sum)[pd_positive_plans]/plan_data[pd_positive_plans,"demand"]
		plan_data[pd_positive_plans,"share_35to54"] <- by(households$weight * (households$perc_35to44 + households$perc_45to54),
			households$unique_id_nocsr,sum)[pd_positive_plans]/plan_data[pd_positive_plans,"demand"]
		plan_data[pd_positive_plans,"share_250to400"] <- by(households$weight * (households$FPL > 2.5 & households$FPL <= 4),
			households$unique_id_nocsr,sum)[pd_positive_plans]/plan_data[pd_positive_plans,"demand"]
		plan_data[pd_positive_plans,"share_gt400"] <- by(households$weight * (households$FPL > 4),
			households$unique_id_nocsr,sum)[pd_positive_plans]/plan_data[pd_positive_plans,"demand"]
		plan_data[pd_positive_plans,"share_male"] <- by(households$weight * households$perc_male,
			households$unique_id_nocsr,sum)[pd_positive_plans]/plan_data[pd_positive_plans,"demand"]
		plan_data[pd_positive_plans,"share_family"] <- by(households$weight * (households$household_size > 1),
			households$unique_id_nocsr,sum)[pd_positive_plans]/plan_data[pd_positive_plans,"demand"]
		plan_data[pd_positive_plans,"share_Asian"] <- by(households$weight * households$perc_asian,
			households$unique_id_nocsr,sum)[pd_positive_plans]/plan_data[pd_positive_plans,"demand"]
		plan_data[pd_positive_plans,"share_Black"] <- by(households$weight * households$perc_black,
			households$unique_id_nocsr,sum)[pd_positive_plans]/plan_data[pd_positive_plans,"demand"]
		plan_data[pd_positive_plans,"share_Hispanic"] <- by(households$weight * households$perc_hispanic,
			households$unique_id_nocsr,sum)[pd_positive_plans]/plan_data[pd_positive_plans,"demand"]
		plan_data[pd_positive_plans,"share_Other_Race"] <- by(households$weight * households$perc_other,
			households$unique_id_nocsr,sum)[pd_positive_plans]/plan_data[pd_positive_plans,"demand"]
		
		plans_pmt$risk_score <- exp(predict(risk_score_reg,newdata=plans_pmt[,fields]))
		plan_data$risk_score <- exp(predict(risk_score_reg,newdata=plan_data[,fields]))

	# Lagged market share
	fields <- c("lagged_market_share","lagged_risk_score","premium_increase_perc","premium_increase")
	plans_pmt[,fields] <- NA
	plan_data[,fields] <- NA
	for(t in 2015:2018) {
		year_plans <- rownames(plans_pmt[plans_pmt$year == t,])
		previously_offered_plans <- paste(plans_pmt[year_plans,"plan_name"],plans_pmt[year_plans,"rating_area"],t-1,sep="")
		
		year_plans <- year_plans[previously_offered_plans %in% rownames(plans_pmt)]
		previously_offered_plans <- previously_offered_plans[previously_offered_plans %in% rownames(plans_pmt)]
		
		plans_pmt[year_plans,"lagged_market_share"] <- plans_pmt[previously_offered_plans,"market_share"]
		plans_pmt[year_plans,"lagged_risk_score"] <- plans_pmt[previously_offered_plans,"risk_score"]
		plans_pmt[year_plans,"premium_increase_perc"] <- 
			plans_pmt[year_plans,"Premium"]/plans_pmt[previously_offered_plans,"Premium"] - 1
		plans_pmt[year_plans,"premium_increase"] <- 
			plans_pmt[year_plans,"Premium"] - plans_pmt[previously_offered_plans,"Premium"]
			
		year_plans <- plan_data[plan_data$Year == t,"unique_id_nocsr"]
		previously_offered_plans <- paste(plan_data[year_plans,"Plan_Name_NOCSR"],plan_data[year_plans,"region"],t-1,sep="")
		
		year_plans <- year_plans[previously_offered_plans %in% rownames(plan_data)]
		previously_offered_plans <- previously_offered_plans[previously_offered_plans %in% rownames(plan_data)]
		
		plan_data[year_plans,"lagged_market_share"] <- plan_data[previously_offered_plans,"market_share"]
		plan_data[year_plans,"lagged_risk_score"] <- plan_data[previously_offered_plans,"risk_score"]
		plan_data[year_plans,"premium_increase_perc"] <- plan_data[year_plans,"Premium"]/plan_data[previously_offered_plans,"Premium"] - 1
		plan_data[year_plans,"premium_increase"] <- plan_data[year_plans,"Premium"] - plan_data[previously_offered_plans,"Premium"]
	}

##### Run regressions

	library(estimatr)
	library(car)

	# Plans pmt
	spec <- premium_increase_perc ~ lagged_market_share + lagged_risk_score + HMO + AV + Anthem + Blue_Shield + Health_Net + Kaiser
	spec <- premium_increase_perc ~ lagged_market_share + risk_score + HMO + AV + Anthem + Blue_Shield + Health_Net + Kaiser
	spec <- premium_increase_perc ~ lagged_market_share + risk_score + HMO + AV + Anthem + Blue_Shield + Health_Net + Kaiser + as.factor(region)
	reg_perc <- lm(spec,data=plans_pmt[plans_pmt$year > 2014,])	
	summary(reg_perc)		
	
	spec <- premium_increase ~ lagged_market_share + risk_score + HMO + AV + as.factor(Insurer)
	reg <- lm(spec,data=plans_pmt[plans_pmt$year > 2014,])	
	summary(reg)	
	
	# Full plan_data (still need to delete CSR plans)
	spec <- premium_increase_perc ~ lagged_market_share + risk_score + HMO + AV + Anthem + Blue_Shield + Health_Net + Kaiser + as.factor(region)
	reg_full <- lm_robust(spec,data=plan_data,se_type = "stata")	
	coef_full <- coef(reg_full)
	se_full <- summary(reg_full)[12][[1]][,2]
	pval_full <- summary(reg_full)[12][[1]][,4]
	summary(reg_full)		
	
	spec <- premium_increase_perc ~ lagged_market_share + risk_score + HMO + AV + Anthem + Blue_Shield + Health_Net + Kaiser
	reg_nomkt <- lm_robust(spec,data=plan_data,se_type = "stata")	
	coef_nomkt <- coef(reg_nomkt)
	se_nomkt <- summary(reg_nomkt)[12][[1]][,2]
	pval_nomkt <- summary(reg_nomkt)[12][[1]][,4]
	summary(reg_nomkt)	
	
	spec <- premium_increase_perc ~ lagged_market_share + risk_score + HMO + AV
	reg_nomktnofe <- lm_robust(spec,data=plan_data,se_type = "stata")	
	coef_nomktnofe <- coef(reg_nomktnofe)
	se_nomktnofe <- summary(reg_nomktnofe)[12][[1]][,2]
	pval_nomktnofe <- summary(reg_nomktnofe)[12][[1]][,4]
	summary(reg_nomktnofe,robust=TRUE)	
	
	
	library(stargazer)
	output_names <- c("Lagged Market Share","Risk Score","HMO","AV")
	stargazer(reg_perc,reg_perc,reg_perc,covariate.labels=output_names,
			style="aer",model.names=FALSE,
			model.numbers=FALSE,column.labels=c("(1)","(2)","(3)"),omit="region",
			coef=list(coef_nomktnofe[2:5],coef_nomkt[2:5],coef_full[2:5]),
			se=list(se_nomktnofe[2:5],se_nomkt[2:5],se_full[2:5]),
			p=list(pval_nomktnofe[2:5],pval_nomkt[2:5],pval_full[2:5]),
			dep.var.labels.include=FALSE,omit.stat=c("f","adj.rsq","res.dev","ser"),
			dep.var.caption="% Change in Premiums")
	